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Abstract 

In this paper we derive accurate numerical methods for the quantum Boltz- 
mann equation for a gas of interacting bosons. The schemes preserve the main 
physical features of the continuous problem, namely conservation of mass and 
energy, the entropy inequality and generalized Bose-Einstein distributions as 
steady states. These properties are essential in order to develop schemes that 
are able to capture the energy concentration behavior of bosons. In addition we 
develop fast algorithms for the numerical evaluation of the resulting quadrature 
formulas which allow the final schemes to be computed only in 0(N 2 log 2 N) 
operations instead of 0(N 3 ). 

Key words: Boson Boltzmann equation, condensation, quadrature formulas, fast 
algorithms. 
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1 Introduction 

We consider a gas of interacting bosons, which are trapped by a confining poten- 
tial V = V(x) with mmV(x) = 0. We denote the total energy of a boson with 
momentum p and position x (after an appropriate non-dimensionalization) by 

I |2 

e( x ,p) = ^- + V(x). (1) 
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Let F = F(p, x,t) > be the phase-space density of bosons. Assuming a boson 
distribution which only depends on the total energy e we write 

F(x,p,t) = f(^ + V(x),t 

where / = f(e,t) > is the boson density in energy space. 

1.1 The Boson Boltzmann equation 

Following [16j,[17j,[18 j ,[19j,[20j we write a Boltzmann-type equation (referred to as 
boson Boltzmann equation in the sequel) in energy space 

p(e)^=Q(f)(e), t>0, (3) 

with the collision integral 

Q(/)(0 = / 6ie + e*-6?-£)S(e,e*,e f ,£)[f , fi(l + f)(l + f,) 
Jr 3 , 



* 5 



(4) 



- + /')(! + fi)]ds,ds'd8 

where S > is a given function. 

We denoted the density of states by 

p(e) = 8 (e - (M! + V ( X )X\ dp dXj ( 5 ) 

and 

f' = f(e',t), fi = f(e[,t), f = f(e,t), /* = /(£*,*). (6) 

As usual e and are the pre-collisional energies of two interacting bosons and e' 
and are the post-collisional ones. 
The positive measure 

5(e + £*-£'- <)£(£, £*,£',<) (7) 

denotes the energy transition rate, i.e. Sde 1 cfe* is the transition probability per unit 
volume and per unit time that two bosons with incoming energies e, e* are scattered 
with outgoing energies e', e*. 

A simple computation shows that the phase-space density F = F(x,p,t) satisfies 
the momentum-position space Boltzmann equation 

OF 

— +p-V x F-V x V(x)-V p F = Q(F), (8) 
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with the scattering integral 

P(N / 2 + K(a;)) 

Note that ([9]) does not correspond to the physical Boltzmann operator for bosons 
except in the homogeneous case V{x) = and F independent of x, where we set 

p(e)= fsL-Wpldp (10) 



and compute 



p{e) = 4irV2e. (11) 



Then equation (|8|) is formally identical to the Boson Boltzmann equation considered 
in [6],[71 



,P',P'*) 



dF f 

— = / 5{p + p* - p ' + £t - e - ejwfap*,! 

[F'Fl{\ + F)(l + F*) - FF*(1 + F')(l + K)\ dp.dp'dp'*, 
with = \p\ 2 /2 and VF, 5 are related by 

S(p + P* - p' - v'*) W (p, P* > P' , P* ) da* da' da* 



(12) 



S 2 xS 2 xS 2 

S(\p\ 2 /2, |^| 2 /2, iyi 2 /2,K| 2 /2) 

p(|p| 2 /2)Ip*IIp'IKI 

Here we denoted p* = \p*\a*, p' = \p'\a', p = \p\a, and p* = Ip'Ja'*. In particular for 
W = 1 we have 

S(e,£*,e',e'*) = const p(e m i n ), (13) 

where (see [6]) 

Smin = min(e, £*,£', 4)- ( 14 ) 

Even in the non-homogeneous case V{x) ^ the equation (|3|) is formally identical 
to the isotropic version of the homogeneous bosonic Boltzmann equation (|12p (after 
the introduction of |p| 2 /2 as new independent variable). However, the density of 
states is computed by formula © in the non homogeneous case instead of (jlOp in 
the space homogeneous case. 

In the physical literature the equation ([3]), usually referred to as ergodic ap- 
proximation of the Boltzmann equation, is derived in the nonhomogeneous case as 
approximation of the phase-space Boltzmann equation by a projection technique 

mm- 
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For a mathematical analysis of the bosonic Boltzmann equation in the space 
homogeneous isotropic case we refer to [10]. [11]. [6]. [7]. We remark that already the 
issue of giving mathematical sense to the collision operator Q(f) is highly nontrivial 
(particularly for scattering rates without cutoff or if positive measures / are allowed, 
as required by a careful analysis of the equilibrium states). 

1.2 Physical properties 

A simple calculation gives the weak form of the collision operator. Let <j> = (j){e) be 
a test function. Then, at least formally 



Here we used the micro-reversibility property, i.e. the fact that each collision is re- 
versible and that each pair of interacting bosons represents a closed physical system. 
Mathematically this amounts to the requirement [16] 



The symmetry properties ([16]) immediately imply the analogous properties for the 
energy transition rate ([7]) and the weak form (|15p follows from the variable substi- 
tution in the integral using these symmetries. 

As a consequence we have the following collision invariants 





S(s + s*-e'- e'*)S(e,E*,e',e:)[f'fUl + /)(1 + /„) 



//*(1 + /')(1 + fl))[<t> + 4>*-4>' ~ <f>i]d£de*d£'dei. 



(15) 



S(£, £*,£',£*) = 5(e*, £,£',£*) = S(e', £*, £,£*). 



(16) 



1. 




(17) 



2. 




(18) 



Consider now the IVP ([3]) supplemented by the initial condition 



/(e,t = 0) = / (e)>0, £ 



> 0. 



(19) 



Then ()17[) implies mass conservation 




(20) 
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and ([18]) energy conservation 

f'OO f'OO 

/ p(e)f(e,t)ede = / p(e)f (e)ede, Vt > 0. (21) 
jo jo 

The H-theorem for ([3]) is derived by setting cp(e) = ln(l + /(e)) — In /(e) in (|15p . 
We calculate 



/■oo 

/ Q(/)(e)(ln(l + /(e))-ln/(e))cfe 
■/ o 

\ f 5(e + e*-e' -£'JS(e,e*,e',£'Mf)dede*de'd£: :=£>[/], 



(22) 



where 
and 



e(/) = *(//.(! + /')(! + fi), /'/:(! + /)(! + /*)) (23) 



z(x, y) = (x — y)(lnx — \ny). (24) 



Since the integrand of the entropy dissipation D[f] is non-negative, we deduce the 
following H-theorem, obtained by multiplying ([3]) by (j)(e) = ln(l + /(e)) — In /(e) 

j t S[f] = D[f], (25) 

which implies that the entropy 

/•oo 

S[f}:= p(e)((l + /)ln(l + /)-/ln/) ( fe, (26) 
J o 

is increasing along trajectories of (|3|). We remark that trivially the third physical 
conservation law, namely momentum conservation, also holds. Clearly the phase- 
space density F of ([2|) satisfies 

/ pF(x,p,t)dx = 0, Vt>0. (27) 

We now turn to the issue of steady states of ([3]). The problem of equilibrium 
distributions for bosons has a very long history, going back to Bose and Einstein in 
the twenties of the last century (see [I] , [1] , [5] ) , who noticed that the class of 'regular' 
Bose-Einstein distributions 

^°( £ ) = e ae+l_^ «>0,/3>0 (28) 
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is not sufficient to assume all arbitrarily large values of equilibrium mass 

roo 

Moo = / p(e)f 00 (e)de, (29) 
Jo 

and arbitrarily small values of equilibrium energy 

;>oo 

£oo= / p{e)ef OQ {e)de, (30) 
Jo 

such that Dirac distributions centered in zero energy have to be included in the set 
of equilibrium states. In [5] it was shown that for every pair (M^, £00) G M? + there 
exist a > 0, /3 £ R such that the generalized Bose-Einstein distribution defined by 

PfcO/eofcO = e J + t}_ l + (31) 

is an equilibrium state of (J3j) (in the sense of maximizing the entropy, see [6] for 
analytical details) satisfying (f29]) - (f30|) . Here we denoted (3 + = max(/3, 0) and /?_ = 
— max(— /3, 0). The value M OC;Con d = |/3_| represents the mass of particles which are 
condensed in equilibrium, i.e. in their quantum mechanical ground state with e = 0. 

Off course, it is analytically nontrivial to define the nonlinearities in the en- 
tropy (dissipation) and in the collision operator, in particular for measures which 
are singular with respect to the Lebesgue measure, as required for the equilibrium 
states. For details we refer to the references [6| and [TO], here we only mention that 
an approximation argument shows that the singular part of a measure / does not 
contribute to the entropy S[f]. For appropriate scattering rates (with unphysical 
cut-off) in the homogeneous case an existence/uniqueness theory for integrable and 
for measure solutions can be set up. So far, it is not clear how the cut-off assumption 
can be removed. 

In the following sections we shall use 

S(e, £*,£',<) = p(e mia ). (32) 

Notice that the condensation is fully localized in phase space, i.e. it may only 
occur at p = (vanishing momentum) and at those points in position space, where 
the potential assumes its minimum value 0. The reason for this is the form ([2]) of the 
phase space distribution and a semiclassical limit process which leads to the Boson 
Boltzmann equation ([3]). 

The purpose of this paper is to derive an accurate discretization of the IVP (|3|), 
(|19h . which maintains the basic analytical and physical features of the continuous 
problem, namely 

• Mass and energy conservation 
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• Entropy growth 

• Generalized Bose-Einstein equilibrium distribution 

To this aim we shall derive first and second order accurate quadrature formulas for 
Q(f). These schemes due to their 'direct' derivation from the continuous operator 
possess all the desired physical properties at a discrete level. In addition we show 
that with the choice (|32[) the computations can be performed with a fast algorithms 
reducing the 0(N 3 ) cubic cost to 0(N 2 log 2 N). For the sake of completeness we 
mention the recent works [2] . [5] . [1^] . |14| . [15] in which fast methods for Boltzmann 
equations were derived using different techniques like multipole methods, multigrid 
methods and spectral methods. 

The rest of the paper is organized as follows. In the next Section we discuss 
the details of our numerical schemes, together with the issues of consistency and 
computational complexity. In Section 3 several numerical tests are performed. The 
results confirm the expected accuracy of the schemes and in particular show the 
ability of the methods to capture the concentration behavior of bosons. Finally we 
concluded the paper with some remarks in Section 4. 



2 Fast, conservative and entropic methods 

We consider the IVP for the quantum boson Boltzmann equation 

p(e)^ = Q(f)(e), t>0, (33) 
/(M = 0) = /oOO>0. (34) 

Here the independent variable e > represents the kinetic energy, p = p(e) > is 
the (given) density of states and the boson collision operator now reads 

Q{f){e) = I 6(e + e*-s'-ei)p(e min )[f'fi(l + + 

(35) 

- //*(! + f')(l + fi)] deeds' de',. 



Obviously the equation (|33p maintains a minimum principle such that solution of 
(1331) . dMD satisfy f(e, t) > for e > 0, t > if / (e) > for e > 0. 



2.1 Reduction on a bounded domain 

Our starting point in the development of a numerical scheme for (|35p is the definition 
of a bounded domain approximation of the collision operator Q. 
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Let / be defined for e E [0, R] and denote 

Q R {f){e) = [ S(e + e,-e/-e'Me n ^)[f'fl(l + + 

J\0,R] 3 



[Q,R] 3 

(36) 

- + /')(! + < R) de.de' de', 

where ip(I) is the indicator function of the set /. Then, at least formally 

l°° QnU)<i>de = -/ 5(e + e*-e> -e>«)p(e min )[ffl(l + + 
2 ^[o,i?] 4 

(37) 

- + /')(! + /*)][«/> + <A* " 0' " 4>'*]dede*de'de'* 

for any test function = ^(e). The proof follows the lines of the corresponding 
weak form of Q discussed in Section 1. 
Consider now the approximate IVP 

p(e)^- = Q R (f R )(e), t>0, (38) 
f R (e,t = 0) = fo, R ( £ )>°- ( 39 ) 
Then the weak formulation (|37p of Q R implies mass and energy conservation 

(40) 



r-R r-R 

/ p(e)f R (e,t)de = / p(e)f . R (e)de, Vi > 0, 
io Jo 



/ p(e)f R (e,t)ede= / p(e)f . R (e)ede, Vt > 0. (41) 
io jo 

Also the entropy inequality 



J t SR[fR] > 0, (42) 



holds with the entropy 



S R [f R ] ■= [ R p(e)((l + f R ) ln(l + f R ) - f R lnf R )de. (43) 
J o 

These properties are in full analogy with the corresponding ones of the IVP ([3]), 

dm 
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2.2 Discretization and main properties 

Let us now introduce the set of discrete energy grid points e\ < £2 < • • • < sn in 
[0, R]. A general quadrature formula for ([36}) is given by 

N 

Q R (f)(ei)^Q R (f)(ei) = Yl W "j^P( £ min)[fkMl + fi)(l + fj) 

j,k,l=l 

(44) 

- / i / j (l + / Jk )(l + /i)]#i<4 

where now /j = /(ei) and e m i n = min{ej, e 3 -, e^, £/}. The quantities wfj are the 
weights of the quadrature formula and Sfj a suitable discretization of the 5-function 
on the grid. 

In order to maintain the conservation properties on the discrete level it is of 
paramount importance that the discretized 5-function will reduce the points in the 
sum to a discrete index set which satisfies the relation i+j = k+l. Thus it is natural 
to restrict to equally spaced grid points which satisfy exactly the aforementioned 
relation on the computational grid. 

We will further simplify the quadrature formula by considering product quadra- 
ture rules with equal weights for which ivfj = WjWkWi = w 3 with w = R/N and 



i=l 



We now consider the set of ODEs which originates from the energy discretization 
of the IVP ||3SD, (139) 

p(e i )^ = Q R (f)(ei), t>0, (45) 
at 

fi(t = 0) = /o,i?(e 4 )>0. (46) 

and prove 

Proposition 2.1 If we define 

§ ki = f l/w i + j = k + l , 47 . 
iJ 1 otherwise 

the solutions of the IVP ft45\ ), ft46\ ) satisfy the following discrete conservation prop- 
erties and entropy principle 

^$>(e*)^fe) = 0, <X £ ) = 1 > 4>{e) =s, (48) 
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N 

E -°(^)^T Z ^ °. Hfi) = (1 + /*) log(l + /;) - /i log (49) 



x . dh(fi 
"' dt 

i=l 



Proof: 



Due to the definition of 5!?J we have the quadrature formula 

N 

Q R (f)(ei) = W 2 p(£min)[fkfl(l + Ml + fj 

3,1=1 
l<k=i+j-l<N 

~ /i/i(l + /fc)(l +//)]• 



(50) 



In particular, for any test function (ft, formula (|50p admits the following discrete 
analogous of the corresponding weak identity for the collision operator 



N , N 



w^QRtfXeiMei) = \w 3 £ p(e mi „)[/ fc /z(l + + f 3 ) 

(51) 



2 



where fa = 4>{£i). The equations (j!5]l are obtained taking 0(e) = 1, and <j){e) = e. 
The discrete entropy inequality can be derived choosing fas) = /i'(/(e)) = ln(l + 
/(e)) — In /(e). In fact, as in the continuous case, we find 

™E^)^- = ^ 3 E pi^ifkMi+mi+fj) 

1=1 i,j,k,l=l 

i+j=h+l 

(52) 

- fifj(l + + flWUi) + h'(fj) ~ h\f k ) - ti(fi)} > 0, 

since 

h'(fi) + h'{fj) - h'(f k ) - h'(fi) = i°g(U + + - iog((i + / fc )(i + fi)W, 

and the function z{x, y) = (x — y)(log x — log y) > for x, y G 
□ 



Remark 1 It is easy to check by direct verification using that these schemes 
admits 'regular' discrete Bose-Einstein equilibrium states of the form 

foofa) = eaei + p _ v ol > 0, p G R. (53) 
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More delicate is the question of 'generalized' discrete Bose-Einstein equilibrium which 
will be discussed later on. 

Remark 2 Clearly one may use other product quadrature rules with different weights. 
However then the definition of a consistent discrete 5-function which satisfies the 
aforementioned conservation laws and entropy principle becomes very difficult. On 
the other hand it is shown in the next section that the choice of quadrature \50\) 
includes numerical methods up to second order accuracy. 



2.3 First and second order methods 

Let us rewrite for e G [0, R] the collision integral (|36p as 

pR pD{e,e') 

QrU){s)= / p(e min )F(£,e',e'«) deeds', (54) 

JO JS(s,e') 

where F( £)£ ',<) = + /)(1 + /*) - + f )(1 + #)], with e, = e' + < - e, 

and S(e,e') = max{e — e' , 0}, D(e,e') = mm{e — e' + R, R} . The integration domain 
for a fixed value of e in the (e',e^) plane is shown in figured) 




Lemma 2.1 We have 



= < 



' Pi?*) 
P«) 



(£',<) el 
(e',e'J ell 

€ IV 



(55) 



where the regions I, II, III, IV represent a partition of the computational domain 
and are shown in figure [IJ 



Region I is characterized by < e' t < e and < e' < e with e' + > e. Thus 
£* = + e' — e < e, £* = £* + £' — £< e', e* = e* + e' — e < e* and hence £ m i n = e*. 

Region II is characterized by e' > e and > e with e' + e£ < R + e. Thus 
£* = e' + — e > e and hence e m i n = e. 

Region III is characterized by R > e' > e and < e' # < e. Thus £* = s'+£* — £ > 
£* and hence e min = £*. 

Region IV is characterized by R > £* > e and < e! < e. Thus £* = £' + £* — £> 
e' and hence £ m i n = e'. 



Using the previous lemma the integral (|54|) over the four regions can be decom- 
posed as 



Proof: 



□ 



QrU){b) = h(e) + h{e) + h{e) + h(e) 



(56) 



with 




(57) 



(58) 



(59) 



(60) 



A similar decomposition holds for the quadrature formula (|50p 



Qn{f){ei) = h(ei) + h{e l ) + IM) + Ufa) 



(61) 
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with 



i i 



life) = w 2 ^2 ^2 p{e k + ei - £i)F(£i,e k ,ei), (62) 

k=l l=i-k+l 
N N+i-k 

I 2 fe) = W 2 ^ ^2 P( £ i) F ( £ i' £ k, £ l)> ( 63 ) 
k=i+l l=i 
N i 

hfe) = w 2 ^ ^2 p(ei)F(£i,e k) ei), (64) 
k=i+i i=\ 

i N 

hfe) = w 2 ^2^2 P{ e k)F{£i,e k ,£i). (65) 
k=i i=i+i 

From the point of view of accuracy we can state 

Theorem 2.1 (Consistency) Let the function f and p be C m ([0, R]), m = 1 or 
m = 2, then the quadrature formula A50\) satisfies 



\Qdf)fe) ~ Q R {f)fe)\ < R 2 C m (Ae) m M m , Ae = R/N, (66) 

where M m is a constant that depends on f and p and their derivatives up to the 
order m and if £% = (i — l)Ae, i = 1,...,N (rectangular rule) then m = 1 and 
C m = 1/2, whereas if £% = (i — 1/2) Ae, i = 1,...,N (midpoint rule) m = 2 and 
C m = 1/24. 

Proof: 

First let us recall the following basic estimate for a composite product quadrature 
rule with equal weights (see [3] for example) 



rb r-d N x N v 

j \ g(x,y)dxdy - AxAy^2^2 g(xi,yj) 

Ja Jc i=l j=\ 



< 



(6 _ a )( d _ c ) Cm [(Ax) m M x , m + (Ay) m M y>m \ , ^ ' 

where Ax = (h - a)/N x , Ay = (d - c)/N y , M x>m and M y , m are two constants such 
that 

d m q, ,d m q, 

y < M < M 

I Q x m I - 1V1 *,rn, \ q | ^ My,m, 

on [a,b] x [c,d] and if X{ = (i — l)Ax, yi = {i—X)Ay then m = 1, C m = 1/2, whereas 
if Xi = (i- 1/2) Ax, yi = (i- 1/2) Ay then m = 2, C m = 1/24. 
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Now, since the integrands which appear in Jj satisfy the required regularity con- 
ditions and approximations given by 1% are the corresponding generalized composite 
product quadrature rules, each error |Jj — Jj| can be estimated similarly to (|67p. 
More precisely we have 

Ih^) - h(Si)\ 

\h{ei) - h{Ei)\ 

\h{Si) - h{Si)\ 
\h{Ei) - h(Ei)\ 

where the constants M l £ , m (e) and M l £ , m (e) are suitable bounds of the partial deriva- 
tives of order m of the integrand functions. 
Summing up the errors we get 

\Q R (f, f)(£i) ~ Q R {si)\ < R 2 C m (Ae) m M m , (68) 

where M m (e) = max^M*, m (e k ) + M\, m (e k )}. 
□ 



< {etfCn^eTWll^Ei) + M £ \ im ( £i )], 

< (R - ^) 2 C m (Ae) m [M e 2 /)m (^) + Mf, tm ( £i )], 

< £i (R - £i )C m (Ae) m [Af £ 3 , Jei) + Ml m (e t )], 

< Ei(R - e i )C m (Ae) m [M £ 4 , m (e t ) + M^Je^ 



2.4 Fast algorithms 

Finally we will analyze the problem of the computational cost of the quadrature 
formula (j50|) . A straightforward analysis shows that the evaluation of the double sum 
in ([50]) at the point £j requires (2(i— 1)(N— i+l)+N 2 ) /2 operations. The overall cost 
for all N points is then approximatively 2iV 3 /3. However using transform techniques 
and the decomposition (|6ip this 0(iV 3 ) cost can be reduced to 0(N 2 log 2 N). 
In order to do this let us set h = k + I = i + j in (j50|) and rewrite 

27V N 



Q R (ei) = v?J2Y,p( e ^fkfh-k(i+Mi+h-i) 

h=2 k=l 



(69) 



where we have set 



*M = J 1 s ^ ^ d (70) 
2 otherwise 



In (j69|) we assume that the function /j is extended to % = 1, . . . , 2iV by padding zeros 
for i > N. 
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The sum (|69p can be split into sum over the four regions which characterize 
p(emin)- We shall give the details of the fast algorithm only for region I, the other 
regions can be treated similarly. We have 

27V i 



h=2 k=l 



(71) 



or equivalently 

27V i 

h=2 k=l 

27V i 

/i=2 fc=l 
27V 

= u, 2 J] p(e fc _i)(l + / f )(l + / A - f )*LM^(i) 

ft=2 
27V 

- w 2 j2p( £ h-i)fif h -Miis 2 h {i), 

h=2 

where we have set 

Sl(i) = £ fkh-k^lX Slit) = (1 + /*)(! + /*-*)*£!■ (72) 



Now the two sums <S^(£) and 5?(i) are discrete convolutions and can be evaluated 
for all /i and i using the FFT algorithm in 0(N 2 log 2 N) operations. This can be 
easily done rewriting them in the form 



TV 

S h (i) = J2^9h-k^i^' i] , (73) 
fc=i 

for a suitable choice of the discrete function gi. It is well known that for N = 2 a 
with a integer the sum (|73|) can be computed for each i via FFT in 0(A r log 2 N)) = 
0(2 a a) operations. The total cost to compute Sh(i) for all i is then 0(N 2 log 2 N). 

A better algorithm can be obtained if we rewrite the sums Sl{i) and S 2 (i) in 
the form 

S h (i) = J29k9 h -k* [ ii^' i] , (74) 
fc=i 
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where 



Pi 



[[\og 2 (i - 1)]] + 1 i > 1 
% = 1 



(75) 



and [[•]] denotes the integer part. 

For each i the convolution sum (|74[) now can be computed in 0(2**/%) operations. 
The total cost will be approximatively reduced by one half since 0(5^i=i *l°g2 ~ 
<9(±iV 2 log 2 (A0). 

Clearly once expressions S^(z) and 5^(i) have been computed the remaining two 
sums are of the type 



2JV 



(76) 



h=2 



which can be computed directly with 0(N 2 ) operations. Thus the final cost for the 
computation of fife) for all i is 0(N 2 log 2 N + N 2 ) = 0(N 2 log 2 N). 



Remark 3 In the case of constant p it is easy to show that expression / fggj) reduces to 
a double convolution sum which can be evaluated using the FFT in only 0(N dog 2 N) 
operations instead of 0(N 2 log 2 iV). 



Figure 2: The relative L\ error for scheme QBF1 (left) and QBF2 (right) computed with 
N = 20 (solid line), N = 40 (dotted line), N = 80 (dash-dot line) points for t € [0, 2.5]. 



3 Numerical tests and applications 

In this section we test the performance of the proposed schemes by considering their 
behavior in different physical and mathematical situations. We shall refer to the 
first and second order fast schemes developed in the previous section by QBF1 and 
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Figure 3: Convergence rates for scheme QBF1 (left) and QBF2 (right) computed with 
N = 20, 40 (solid line), N = 40, 80 (dotted line) points for t £ [0, 2.5]. 



QBF2 respectively. The time integration is performed with standard first and second 
order explicit Runge-Kutta schemes after dividing equation ([43]) by p{Ei) and thus 
rewriting the semidiscrete schemes as 



l<k=i+j-l<N 



In all our numerical tests the density of states is given by 



(77) 



P(e) = y, (78) 

which corresponds to an harmonic potential V(x). 

Note that < jo(e m in) / p( £ i) < 1 for £j / and that as £j — > we have 
p(emin) / 'p( £ i) ~^ 1- Furthermore since p(0) = the values of the distribution func- 
tion at £j = does not affect the discrete conservation of mass and energy. 

The schemes were implemented using the fast algorithm described in Section T2.4I 



3.1 Accuracy analysis 

The first test case has been used to check the numerical convergence of our quadra- 
ture formulas by neglecting the time discretization error (as usual this can be 
achieved either using very small time steps or sufficiently accurate time discretiza- 
tions). The initial datum is a Gaussian profile centered at R/2 

/ = exp(-4(e- J R/2) 2 ), (79) 
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with R = 10. The final integration time is T = 2.5. We report in Figure[2]the relative 
errors in the Li— norm obtained with the different schemes for ./V = 20,40,80 grid 
points. As a reference solution we used the numerical result obtained with a fine 
grid of N = 160 points. 

In Figure [3] the corresponding convergence rates of the schemes are reported. As 
usual given two error curves En and E^iv corresponding to N and 2N grid points 
the convergence rate is computed as 




The results confirm the expected first order and second order degree of accuracy of 
the methods. 

Remark 4 Since the midpoint rule, similarly to the trapezoidal rule, admits an 
Euler-MacLaurin expansion we can in principle increase the order of the method 
by extrapolation techniques. Unfortunately with this approach it is difficult to keep 
conservations as well as entropy inequality. 

3.2 Bose-Einstein equilibrium 

Next we consider the same initial data as in the previous section and compute the 
large time behavior of the schemes for N = 40. The stationary solution at t = 10 is 
given in Figure 0] for both schemes together with the numerically computed entropy 
growth. As observed the methods converge to the same stationary state given by a 
'regular' discrete Bose-Einstein distribution. 




0123456789 10 0123456789 10 



Figure 4: Stationary discrete Bose-Einstein equilibrium and entropy growth for scheme 
QBF1 (o) and QBF2 (x) computed with N = 40 points. 
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The trend to equilibrium in time for the two schemes is reported in Figures [5j 
Note that although the two schemes agree very well there is a remarkable resolution 
difference in proximity of the point e = due to the staggered grids of the schemes. 





Figure 5: Trend to equilibrium in time for scheme QBF1 (left) and QBF2 (right) computed 
with N = 40 points. 

However since the value of / at e = does not affect the macroscopic quantities 
and the entropy we can adopt a suitable extrapolation strategy to recover a better 
resolution of scheme QBF2 near e = 0. Since we are mostly interested in the large 
time behavior of the solution we can recover the value at the zero energy level by 
a steady state extrapolation. This corresponds to assume / of the form (|53|) and 
consequently to assign 

In Table 13.21 we compare the extrapolated results at the final computation time 
of scheme QBF2 for different extrapolation methods with scheme QBF1 and with 
the "exact" steady state solution. We remark that the values of a and f3 for the 
stationary state can be computed by inverting numerically the equations (|29p - (j30p 
for /oo given by (|28|). The marked improvement in the resolution given by scheme 
QBF2 with steady state extrapolation is evident. 

In Figure [6] we present the corresponding result for scheme QBF2 with steady 
state extrapolation at e = (as we shall always do from now on with QBF2). In the 
same figure we also report the final "steady" solution at t = 10 for the phase-space 
density reconstructed at x = and p = (pi,p2,0). 
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Exact 


QBF1 


QBF2 with extrapolation 






Steady state 


Exponential 


Cubic 


Linear 


7.144 


6.335 


7.217 


6.449 


6.323 


5.994 



Table 1: Values of /(0) at t = 10 with TV = 40 points. 




Figure 6: Trend to equilibrium in time for scheme QBF2 (left) and stationary phase-space 
density reconstructed at x — and p — (pi,P2jO) (right) with steady state extrapolation at 
e = 0.' 



3.3 Condensation 

In this test we consider the process of condensation of bosons. It is a fundamental 
results of quantum statistics of bosons that above a critical density/below a critical 
energy particles enter the ground state, i.e. a Bose-Einstein condensate forms (see 
[19], [20], [18], [16], [IT]) and the equilibrium distribution is of the form (f3Tj) with 

In general the evaluation of the condensate fraction as a function of time is a 
challenging problem from the computational viewpoint. If we assume the density 
function / to be of the form pip , which corresponds to the long time behavior, we 
can use the following method to identify if condensation will occur and compute the 
equilibrium condensate mass for a given mass energy pair (M, E) . 

First solve numerically for a the equation 

E.f°-!^-^. (8!) 
J exp(ae) - 1 
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Then compute 

-ds. (82) 



/ 

Jo 



exp(ae) — 1 

If 1(a) < M the mass entropy pair is critical and condensation will take place. The 
condensate mass fraction in equilibrium can then be computed 

M c I a 

We report in Figure [7] the condensate mass fraction computed with the previous 
method for (M, E) e [0, 1] x [0, 1]. 




Figure 7: Mass fraction of the condensate in the mass-energy plane at the stationary state. 

A related challenging problem is the computation of the critical time at which 
the condensate starts to form. In order to do this we consider two different numerical 
indicators. 

We recall that for the second order method, unlike the first order one, due to the 
midpoint quadrature, we have e% ^ for all gridpoints. This makes scheme QBF2 
more suitable to treat situations where the solution is close to be singular at e = 0. 
In particular, in such cases, it is impossible to extrapolate the value /(0) with a 
positive [3. Thus whenever steady state extrapolation is impossible we can assume 
to have formation of condensate at e = 0. 

For the scheme QBF1 we expect the value of /(0) to increase dramatically when 
formation of condensate takes place. In this case we can use as an indicator of the 
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formation of condensate the expression |13j 



C F 



(84) 



For the numerical test we choose the initial distribution in the energy interval 
[0, R] with R = 10 to bepl],[17] 



/(e) = — arctan(e 

7T 



r(l-e/eo)> 



(85) 



with r = 5 and eo = R/8. At values of / larger than a critical /* the formation of 
a condensate occurs (see [E],[T7] for similar results in the homogeneous case). We 
choose / = 1, which turns on to be supercritical. In this case the mass energy pair 
is approximatively (0.42, 0.50) which corresponds to a condensate mass fraction of 
« 0.3 at the stationary state (see Figure [7]). Using N = 320 points and scheme 
QBF2 with steady state extrapolation the condensate formation in finite time at 
t r ~ 4.2 is observed. 



- N=320 
N=160 
N=S0 





Figure 8: Estimation of the critical time using the numerical indicator Cf in time for 
N = 80, 160, 320 for scheme QBF1 (left) and scheme QBF2 with steady state extrapolation 
(right). 

We report in Figure [8] the time evolution of the indicator (|84p for scheme QBF1 
and for scheme QBF2 with steady state extrapolation before the critical time. The 
vertical line correspond to the critical time at which the steady state extrapolation 
fails. The results indicate the numerical convergence of the approximation (|84p . 

The distribution of bosons at different times in logarithmic scale before the criti- 
cal time is shown in Figure [9] for scheme QBF1 (left) and scheme QBF2 with steady 
state extrapolation (right) . 



22 



Figure 9: Distribution of bosons at different times in logarithmic scale before the critical 
time for scheme QBF1 (left) and scheme QBF2 with steady state extrapolation (right) for 
N = 320. 

A magnified view of the numerical solutions obtained with N = 40 and N = 80 
points at t = 15 shows that away from the singularity the two schemes are still in 
good agreement (see Figure [T0|) , 

Finally in Figure [TT] we also report the phase-space density reconstructed at 
x = and p = (pi,P2,0) at two different times before the critical time. The cor- 
responding solution has been obtained for N = 80 with scheme QBF2 and steady 
state extrapolation. 

4 Conclusions 

We have developed first and second order fast solvers for the Boson Boltzmann 
equation assuming a boson distribution which only depends on the total energy. 
The methods preserve all the relevant physical properties (conservation of mass and 
energy, entropy inequality and steady states). The performance of the schemes 
has been tested for both Bose-Einstein and generalized Bose-Einstein steady states. 
The numerical methods have shown the capability to describe well the challenging 
phenomenon of condensation of bosons. 

We remark that, to our knowledge, this is the first example of accurate, conser- 
vative and fast deterministic numerical method for a Boltzmann equation. Previous 
results were available in the literature for Fokker-Planck-Landau type equations (see 
]2],[8],[14J) or using some suitable approximations of the Boltzmann equation (sec 
the recent review [12] and the references therein). 

Note that the present numerical methods can be applied directly even to the 
case of the energy dependent quantum Boltzmann equation for Fermions as well as 
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Figure 10: Magnified view of the distribution of bosons after the critical time of conden- 
sation with scheme QBF1 (o) and scheme QBF2 (x) with N = 40 (left) and N = 80 (right) 
points at time t = 15. 

the classical Boltzmann equation of rarefied gas dynamics. 

We hope to extend in the future these ideas to time dependent potentials [T3] . 
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